/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include "estep.h"
#include "prating.c"

static void
safe_exponentiate_normalize (double*      x,
                             unsigned int n)
{
  /*
   * ok we want to compute 
   * 
   *            exp (x[i])          exp (x[i] - max)
   * x[i] = ----------------- = ------------------------
   *        sum_i exp (x[i])     sum_i exp (x[i] - max)
   *
   * the RHS is more numerically stable
   */

  unsigned int maxindex = 0;

  for (unsigned int i = 1; i < n; ++i)
    {
      if (x[n] > x[maxindex])
        {
          maxindex = n;
        }
    }

  double denominator = 0;
  double max = x[maxindex];

  for (unsigned int i = 0; i < n; ++i)
    {
      x[i] = exp (x[i] - max);
      denominator += x[i];
    }

  for (unsigned int i = 0; i < n; ++i)
    {
      x[i] /= denominator;
    }
}

void
estep (const double*            alpha,
       unsigned int             n_labels,
       const double*            logbeta,
       unsigned int             n_images,
       double*                  pz,
       const double*            priorz,
       const Rating*            ratings,
       unsigned int             n_ratings,
       int                      symmetric)
{
  double logpriorz[n_labels];

  for (unsigned int i = 0; i < n_labels; ++i)
    {
      logpriorz[i] = log (priorz[i]);
    }

  for (unsigned int i = 0; i < n_images; ++i)
    {
      memcpy (pz + i * n_labels, logpriorz, n_labels * sizeof (pz[0]));
    }

  for (unsigned int i = 0; i < n_ratings; ++i)
    {
      for (unsigned int j = 0; j < n_labels; ++j)
        {
          auto double alphafunc (unsigned int);
          double alphafunc (unsigned int k)
            {
              if (k == j) { return 0.0; }
              else if (symmetric && k < j) 
                {
                  return alpha[ratings[i].worker * n_labels * n_labels +
                               k * n_labels +
                               j]; 
                }
              else 
                { 
                  return alpha[ratings[i].worker * n_labels * n_labels +
                               j * n_labels +
                               k]; 
                }
            }

          pz[j + ratings[i].image * n_labels] += 
            logprating (ratings[i].rating,
                        alphafunc,
                        logbeta[ratings[i].image],
                        n_labels);
        }
    }

  for (unsigned int i = 0; i < n_images; ++i)
    {
      safe_exponentiate_normalize (pz + i * n_labels, n_labels);
    }
}
